%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Instable Fermi-Bose mixture [Chui04].

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

Constant;

col='krgbmcy';

figure(1);					clf;

txt1='set(gca,''xLim'',10.^[0 1],''yLim'',[0 1]*1e3,''LineWidth'',1.4,''FontSize'',16);';

txt2='xlabel(''\fontsize{16}{\it\sigma}   (\mum)'');';

txt3='ylabel(''\fontsize{16}{\it\mu_\alpha}   (h_P)'');hold on;';



sgm=10.^[-7:.01:-2];		alf=1/2./sgm.^2;								%condensate size

N  =[1:.2:2]*1e3;			wR = 2*pi*50;									%atom number & trap frequency

m  =    7*uA;																%mass

as =-27.3*aB;				g  =4*pi*hP^2*as/m;								%boson-boson scattering length

r  =[.01:.01:1]*50e-6;		U  =m/2*wR^2*r.^2;								%radial coordinate & harmonic potential



for jj=1:length(N)

% variational method.

	HA =3*hP^2*alf/2/m+m*wR^2/2./alf+N(jj)*g.*(alf/pi).^(3/2);



% numerical integration.

	for j=1:length(alf)

		phi =(2*alf(j)/pi)^(3/4)*r.*exp(-alf(j)*r.^2);						%wave function

		Heff=-hP^2/2/m*(-6*alf(j)+4*alf(j)^2*r.^2).*phi.^2 + U.*phi.^2 + N(jj)*g*phi.^4./r.^2;	

		HB(j)=4*pi*sum(Heff)*(r(2)-r(1));

	end

	subplot(221),semilogx(sgm*1e6,HA/hP,col(jj),'LineWidth',1.4);	eval([txt1,txt2,txt3]);

	subplot(222),semilogx(sgm*1e6,HB/hP,col(jj),'LineWidth',1.4);	eval([txt1,txt2,txt3]);

end

